function [V,D] = jacobi1(A,epsilon,show)
% [V,D] = jacobi1(A,epsilon,show)
% To compute the eigenpairs of a symmetric matrix.
% The original Jacobi`s method of iteration is employed.
% A is an n x n matrix, input.
% epsilon is the tolerance, input.
% V is the diagonal matrix of eigenvectors, output.
% D is the diagonal matrix of eigenvalues,  output.
if nargin==2, show = 0; end
D = A;
[n,n] = size(A);
V = eye(n);
done = 0;
working = 1;
stat = working;
cntr = 0;
[m1 p] = max(abs(D-diag(diag(D)))); % Select element
[m2 q] = max(m1);                   % element of largest
p = p(q);                           % magnitude.
while (stat==working),
  t = D(p,q)/(D(q,q) - D(p,p));
  c = 1/sqrt(t*t+1);
  s = c*t;
  R = [c s; -s c];
  D([p q],:) = R'*D([p q],:);
  D(:,[p q]) = D(:,[p q])*R;
  V(:,[p q]) = V(:,[p q])*R;
  cntr = cntr+1;
  if show==1,
    home; if cntr==1, clc; end; 
    disp(['Jacobi iteration No. ',int2str(cntr)]),disp(''),...
    disp(['Zeroed out the element  D(',num2str(p),',',num2str(q),') = ']),...
  	 disp(D(p,q)),disp('New transformed matrix  D = '),disp(D)
  end
  [m1 p] = max(abs(D-diag(diag(D))));
  [m2 q] = max(m1);
  p = p(q);
  if (abs(D(p,q))<epsilon*sqrt(sum(diag(D).^2)/n)), stat = done; end
end
D = diag(diag(D));
